%% ==================== [Find the statistic of the estimated output] ========================
clear
close all
polOrd=6; %Order       
inVarMat=n1PolyRoots(polOrd); % matching point matrix, n+1th order polynomial root plant as matching point

outVarMat=112+3*(inVarMat*2).^2+exp(inVarMat)-0.2*inVarMat.^3 ...
+0.5*sin(inVarMat);%original system expression
samplePointNum=200000; %sampling points
                                                                     
[PCmean,PCvar,MC_PCmean,MC_PCvar,MC_PCdata,MCinput]=GetStatistics ...
    (inVarMat,outVarMat,polOrd,samplePointNum);
syms x
 
MCoutput=112+3*(x*2).^2+exp(x)-0.2*x^3+0.5*sin(x);
MCoutput=subs(MCoutput,sym('x'),MCinput); 
MCoutput=double(vpa(MCoutput,4));
[Xpdf,X] = ksdensity(MC_PCdata); 
MCmean=mean(MCoutput);
MCvar=var(MCoutput);
[X1pdf,X1] = ksdensity(MCoutput); % find the probability density function of the random variable
%% ========================== [plot] =================================
plot(X1,X1pdf,'LineWidth',1)
xlim([100,200])
hold on
plot(X,Xpdf,'---','LineWidth',1);
title('Probability density of random variables')
xlabel('random variable Y')
ylabel('probability density \it{f}')
legend('MCpdf','PCpdf')




    